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The dynamics of M-site, iV-particle Bose-Hubbard systems is described in quantum phase space 
constructed in terms of generalized SU(M) coherent states. These states have a special significance 
for these systems as they describe fully condensed states. Based on the differential algebra developed 
by Gilmore, we derive an explicit evolution equation for the (generalized) Husimi-(Q)- and Glauber- 
Sudarshan-(P)-distributions. Most remarkably, these evolution equations turn out to be second 
order differential equations where the second order terms scale as 1/N with the particle number. 
For large N the evolution reduces to a (classical) Liouvillian dynamics. The phase space approach 
thus provides a distinguished instrument to explore the mean-field many-particle crossover. In 
addition, the thermodynamic Bloch equation is analyzed using similar techniques. 

PACS numbers: 03.75.Lm, 03.65.-w 
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The phase space formulation of quantum mechanics is 
nearly as old as the theory itself [l|. Although the rep- 
resentation is equivalent to the Schrodinger or Heisen- 
berg picture the resemblance between the classical and 
the quantum phase space description reveals interesting 
analogies and differences between the two regimes. 
However, the usefulness of this approach is by no means 
restricted to illustrations. In quantum optics there is a 
wide range of applications of phase space methods (for a 
general overview see, e.g., 0). One particular technique 
which we will exploit in this paper is the association of 
non-commuting operator equations with c-number differ- 
ential equations. In the case of the position or momen- 
tum representation this differential form of the operators 
is common knowledge. By the same token the correspon- 
dence between an operator acting on a density operator 
and a differential operator acting on a phase space dis- 
tribution in flat phase space is widely used, e.g. in the 
context of quantum noise [3[. Strangely enough, these 
methods were for a long time restricted to the description 
of systems which can be described by the dynamic group 
of the harmonic oscillator like a spinless non-relativistic 
quantum particle or a mode of the quantized radiation 
field. Only eight years ago a general algorithm to con- 
struct an s-parametrized family of phase space distribu- 
tions for systems with arbitrary dynamical Lie groups has 
been proposed Q ■ Therefore it has taken thirty years to 
extend the work of Cahill and Glauber [f| and Agarwal 
and Wolf Q from the Heisenberg-Weyl group to phase 
space topologies differing from the complex plane. 
In this paper, we will present a phase space analysis of 
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with dj , dj being the bosonic annihilation and creation 
operators. This model is a paradigm for the study of 
strongly correlated bosonic systems, describing two ap- 
parently very different systems, Josephson junction ar- 
rays and bosons in optical lattices (see, e.g., [7J and refer- 
ences therein). In both cases, the parameter U describes 
the on-site interaction between the bosons, the hopping 
element A gives the tunneling strength confined to near- 
est neighbors, and tj represents the chemical potential 
at each site j. In dependence of the parameter ratio, the 
system undergoes a quantum phase transition from a su- 
perfluid phase for A ^> U, characterized by long range 
coherence and vanishing gap in the excitation spectrum, 
to the Mott phase for U ^> A, dominated by localization 
effects (§}. Especially the prediction Q and the spec- 
tacular experimental realization [10( of the latter system 
attracted a lot of interest, since this shows that optical 
lattices can be seen as a kind of laboratory for strongly 
correlated many-body systems. 

The dynamical group of the Bose-Hubbard model for M 
sites is spanned by the normally ordered operators djdfc 
with j, k E {1,2, M} and is hence equivalent to the 
special unitary group SU(M). This is underlined by the 
fact that every group element as well as the Hamilto- 
nian itself commutes with the particle number operator 
N = Ylj=x OjCij. Consequently an analysis in terms of 
the flat phase space and the use of related methods, like 
Glauber coherent states, is not adequate. For instance, 
the single operators a,j , dj lead to Hilbert spaces with 
different particle numbers and the order parameter (dj) 
obviously vanishes. These facts have been taken into ac- 
count by some recent approaches [111, [HJ . 
In the present paper we will show that taking into con- 
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sideration the particle number conservation explicitly has 
significant advantages regarding the physical interpre- 
tation and the justification of common approximations: 
Since the dynamical group is no longer a direct sum of 
the Heisenberg-Weyl group, but given by SU(M) sym- 
metries, one has to apply an extended concept of coherent 
states [13( | . These states obey a generalized minimum 
uncertainty relation and stay coherent under an evolu- 
tion which is linear in the generators of the dynamical 
group. Moreover, as we will argue in this paper, the 
corresponding generalized coherent states are equivalent 
to the fully condensed states and are therefore of high 
physical significance. Thus, an analysis in terms of phase 
space distributions based on these states emphasizes di- 
rectly every deviation from a product state matching a 
macroscopic wave function. These states are the basis 
of the approximate description by the discrete Gross- 
Pitaevskii equation, which qualifies the phase space dis- 
tributions as an excellent tool to analyze and illustrate 
the mean-field many-particle correspondence. Further- 
more, the presented method conserving the SU (M) sym- 
metry is particularly suitable to derive and justify mean- 
field equations and truncated phase space approaches. 

The paper is organised as follows: In the next sec- 
tion we will recapitulate the concept of generalized co- 
herent states and discuss the relevant cases. Here we 
will also show that every condensed states can be writ- 
ten as a SU(M) coherent state and vice versa. In the 
third section, we will introduce a method to map opera- 
tor equations onto c-number differential equations for the 
SU(M) algebra. This technique will enable us to calcu- 
late the exact phase space dynamics for the Husimi-(Q)- 
and the Glauber-Sudarshan-(P)-function of the Bose- 
Hubbard model which is without any approximations or 
restrictions to the initial state given by a second order 
linear differential equation in the parameter space of the 
SU(M) coherent states. A comparison to the classical 
Liouville equation in phase space reveals a deeper connec- 
tion: The exact phase space dynamics consists of a first 
order differential equation plus a many-particle quan- 
tum correction of second order decaying with the particle 
number as 1/N. This yields an obvious justification for 
a truncation of the evolution equations for large parti- 
cle numbers, in contrast to established methods as the 
truncated Wigner approach [3], where the justification 
is rather difficult. The first order differential terms can 
be thought of as a classical term since they are identi- 
cal to the results of the Liouville equation. However, this 
technique is not restricted to dynamics. As another possi- 
ble application, we will map the thermodynamical Bloch 
equation onto a differential equation. Finite tempera- 
ture effects in the Bose Hubbard model as, e.g., thermal 
fluctuations have recently attracted a lot of experimen- 
tal and theoretical interest [H, [ijj • A closer analysis 
shows that the resulting density matrix can be also de- 
composed into a classical contribution, affected only by 
the Gross-Pitaevskii Hamiltonian function, plus a many- 
particle correction. These examples show that the phase 



space approach is a distinguished instrument to explore 
the mean-field many-particle crossover. 

II. GENERALIZED COHERENT STATES 

The basic ingredient which we will need in the 
following is the concept of generalized coherent states 
for systems with an arbitrary dynamical Lie group [T~3T ] . 
The parameter space of the generalized coherent states 
determines the corresponding phase space and reflects 
the physical properties of the system by its geometric 
structure. Moreover it has been shown that one can 
construct explicitly a family of phase-space distributions 
for a system with arbitrary Lie group symmetry relaying 
on this concept Q. In this section we will provide the 
basics and the notations for the following. 

So, let G be the dynamical Lie group of the relevant 
quantum system. For simplicity we assume that G is 
connected, simply connected and has a finite dimension, 
which is the case for the matrix Lie groups considered 
in this paper. It is important to note that the general 
approach does not rely on these assumptions. The uni- 
tary irreducible representation of the dynamical group G 
acting on the Hilbert space will be denoted by T. With 
these preliminaries, we can define the generalized coher- 
ent states by the action of an element of the unitary irre- 
ducible representation T on a fixed normalized reference 
state |^o) : 

|^>=T(0)hto>, geG. (2) 

Even though the choice of the reference state is in 
principle arbitrary, it influences strongly the shape of the 
coherent states and the structure of the corresponding 
phase space. Therefore a physically motivated choice 
would be an extremal state of the Hilbert space like 
the vacuum ground state for the Heisenberg-Weyl 
group or the lowest/highest spin state for the case of 
SU(M). Mathematically these states correspond to the 
highest/lowest weight states of the unitary irreducible 
representation [l8l ]. 

The isotropy subgroup or maximum stability group 
H C G consists of every element which leaves the ref- 
erence state invariant up to a phase factor. Formally one 
can write 

T(h)\ip ) = e^l^o) with <j)(h) el VheH. (3) 

With respect to the coherent states, there is a unique 
decomposition for every element g € G into a product of 
two elements, one of the isotropy subgroup H and one of 
the coset space G/H: 

g = Qh, geG, he H and Q e G/H. (4) 

Hence, there is a one-to-one correspondence between the 
elements 51(g) of the coset space H/G and the coher- 
ent states |f2) = \ipn) which preserves the algebraic and 
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topological properties. This construction guarantees the 
characteristic property of the coherent states: a coherent 
state stays coherent under a time evolution linear in the 
generators of the dynamical group. 

Another important property we will need in the follow- 
ing is the (over) completeness of the coherent states fl3j . 
which leads to the resolution of the identity operator of 
the Hilbert space, 



|n> (n|dy*(n) =i, 



(5) 



G/H 



where d/i(f2) denotes the invariant measure on the coset 
space. Moreover, this fact guarantees that one can 
uniquely reconstruct the density matrix from the the P— 
or Q-distribution. 



A. Glauber states 

As the Wigner function [l[ and the Moyal quantization 
[l9j |. the coherent states where first introduced for the 
Heisenberg-Weyl algebra hi = {a, a' ,a'a = n,I}, with a 
and a' being the bosonic annihilation and creation oper- 
ators. One of the first applications was the description of 
a mode of the quantized radiation field modeled by har- 
monic oscillators [20| . In this case the unitary irreducible 
representation of an arbitrary group element g £ can 
be decomposed as 



T(g) 



a E C, <5, £ 



(6) 



with the stability subgroup U(l) x U(l) being generated 
by {n, I}. Therefore the phase space is isomorphic to the 
complex plane Hi/U{l) x U(l) = C, parametrized by 
the complex parameter a and the typical representative 
of the coset space 



D(a) = 



(7) 



is just the well-known displacement operator. With the 
physically motivated choice of the vacuum ground state 
|0) as the reference state one obtains the famous Glauber 
states 



D(a)\0). 



(8) 



The generalization to more then one mode is straight- 
forward, since the multimode group © ieN {dj, a\, ajai = 
hi, 1} is just a direct sum of the single-mode group. Thus 
the multimode Glauber states can be obtained as a direct 
product of the single-mode Glauber states, 

M 



n 
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M 
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with |0) being the multimode vacuum ground state. Due 
to this factorization the well-known properties of the 
single-mode Glauber states can be transferred easily. 



B. S£/(M)-coherent states 

In the case of the Bose-Hubbard model {1} with M 
sites, the dynamical group is equivalent to the special 
unitary group SU(M), spanned by the generalized an- 
gular momentum operators Ejk = ajdfc with j, k £ 
{1,2,..., M}. These fulfill the algebraic commutation re- 
lations 



j E mn 



— Ejn&km E ra \ t b 
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and conserve the particle number N = _ 1 Ejj, since 



(10) 



E jk ,N =0 



(11) 



As already argued above, a suitable choice of the ref- 
erence state is the maximum spin state, corresponding 
to the state with the entire population in the first well 
\N, 0, . . . , 0). With respect to this state, an arbitrary ele- 
ment of the unitary irreducible representation can always 
be decomposed as 

T(g)\N, 0, . . . , 0) = exp {^{y k iE kl + yi k E lk )j 
x exp | VkiEu + VixEu ) |JV, 0, . . . , 0) (12) 

\k,l=2 J 

into an element of the coset space and an element of 
the stability group U(M - 1) x U(l) [Hj]. Given that 
Ejk = E^p we have to assume that y* k = y k j in order 
for the argument of the exponentials to be anti-hermitian. 
Therefore we get the SU(M) coherent states by the ac- 
tion of the representative of the coset space onto the ref- 
erence state 

£(y)|jV,0,...,0) 



M 



exp (J2(ykA-yki E lt) |iV,0,...,0) 



k=2 



y)- 



(13) 



The parameter space of the coherent states is spanned 
by the M — 1 complex parameters yk = j/fci with k £ 
{2, ...,M} of the coset space and can thus be identified 
with the 2(M — 1) sphere which is topologically equiv- 
alent to U(M)/U(M - 1) x [7(1) = SU(M)/U{M - 1). 
Due to this analogy one can interpret the coset represen- 
tative as a rotation of the reference state on the multi- 
dimensional sphere. To assure that the parametrization 
is unique one has to demand that the parameters are 
bounded as J2h=2 DkVk < (^/ty 2 - In the case of two sites 
the definition of the coherent states reduces to the spin 
coherent states or Bloch states [H, H3 • 

Anyhow, a parametrization by the [M — 1) indepen- 
dent complex parameters (x2,---,Xm) of the site to- 
gether with the real dependent parameter of the first site 
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x* = x\ is physically more reasonable. These parame- 
ters represent the probability amplitudes at the respec- 
tive sites, reflect directly the particle conservation 



M 

X l + ^A X 2 = 1, 

k>2 



(14) 



and the irrelevance of the global phase. By means of the 
generalized Baker-Campell-Hausdorff formula one can 
show the relation 



Kalll- 1 = cos(||j/||)oJ + 



-t , sin (bll) 



A I 



V 



(15) 



k=2 



with the abbreviation \\y\\ z = Ylk=2 \Vk\' 2 - This leads 
directly to the parameter transformation 

Xl =008(112/11), x k = ^Mly k , k>2 (16) 



I. VI 



and the representation of the SU{M) coherent states in 
terms of the complex amplitudes (xi, x%, . . . , Xm)' 

|y) = n\N,0,...,0) 



fNl 



7ear|o,o,...,o) 



1 / M \» 

i ( M \ N 

7^(X>4) |o,o,...,o) 



X N, 



(17) 



where we have used the commutation relation (|15p . The 
last relation reveals another interesting property of the 
SU(M) coherent states. In the case of the Bose-Hubbard 
model these states are equivalent to the fully condensed 
states, since they can always be written as a product 
state. This characteristic trait is certainly not trivial and 
it cannot be generalized to other dynamical groups since 
it is an intrinsic property of the su{M) algebra. More- 
over, this fact also singles out the physical significance 
of an analysis in terms of phase space distribution which 
are based on the SU (M) coherent states. 



III. DIFFERENTIAL ALGEBRA 

In this section we will present a formalism to map 
quantum observables onto differential equations acting 
on the continuous parameter space of the coherent states 
based on the ideas of Gilmore [2l| which we will use 
to calculate the exact phase space dynamics for the 
Bose-Hubbard model. In contrast to other approaches, 
for example based on the star product (see [24| and 
references therein), this formalism is not restricted to 
the case of just two sites or to the special case of some 
dynamical groups [l8j . 



Flatland 



In the field of quantum optics the modus operandi for 
the Heisenberg-Weyl group H4 and the Glauber coherent 
states is well-known (see, e.g., Q and references therein). 
Since the Glauber states expressed in Fock states \n), 



n— 



n=0 

In' 



(18) 



form an overcomplctc basis, one can replace the action of 
the bosonic creation and annihilation operators by first 
order linear differential equations acting on the function 
f n (a) = exp(— \aa*)u n j\pn\. This yields the differential 
operators V k acting on a ket state 



A\a) = V k {A)\a) 

with V k (a t ) = ^- + la* and V k (a) 
aa 2 



(19) 



Since we are in the following interested in phase space 
densities corresponding to density operators and there- 
fore to products of functions f n (a)f m (a*), we need the 
differential operators T> 1 acting from the left side on the 
coherent state projectors: 



A\a) (a\ = V l {A)\a) (a\ 







with £>'(a T ) = — 
da 



and V l (a) = a. (20) 



The generalization to operators V acting from the right, 



v r {A) = 



(21) 



and to multimode Glauber states is straightforward: 

OOLi 



(22) 



By means of the properties of the differential operators 
acting on arbitrary elements of the multimode algebra 
A, B with r, s e C, 



V l {rA + sB) = rV l {A) + sV l (B) 
V l {AB)=V l {B)V l {A) 



A, B 



V l (B),V l {A) 



(23) 
(24) 

(25) 



one can show that the differential operators conserve the 
algebraic structure. Therefore the differential operators 
of the generators of the Heisenberg-Weyl algebra form 
itself a closed (differential) algebra. 
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B. From the plane to the sphere 



A short calculation gives the desired result 



The su{M) algebra is generated by the set of oper- 
ators {Ejk = a]afc} with j, k e {1, 2, 3, M}. In the 
case of the multimode Glauber states, the corresponding 
differential operators read 



V l (E jk ) = V l {a k )V\a]) = a k d aj + a k a* 
Using the transformation 



XiCte 



ail 



(26) 



(27) 



to the M — 1 complex parameters x = (x2, X3, ■ ■ • , %)', 
the norm a and the global phase <j>, one obtains the dif- 
ferential form of the generalized angular momentum op- 
erator in terms of the multimode Glauber states 



d 



a d 



V l {E jk ) = x k —+x k x- l -- 

dxj 3 \ 2 da 

-ix fc a$(xV + x*V*). 



(28) 



Here we have used the definition 

M 



xV + x*V 



E 



d 
dxk 



_d_ 

'dxt 



(29) 



k=2 " k 

The parameter x\ = x\ is fixed by the normalization 



Xi 



\ 



M 

E x k x ki 

k=2 



(30) 



= X> z (4- fc )|x) J v(x| JV , (35) 
where we used the following abbreviation: 

d 1 , . 9 , , 

^-^ xV - xV ^-^- (36) 

A comparison to equation (f3"Tj) shows that the differenti- 
ation no longer depends on the global phase. This can be 
understood as an averaging effect of the integration over 
the angle 0, which is part of the homomorphism. 

IV. DYNAMICS 

A. The Husimi— distribution 

The time evolution of the Husimi- or Q-distribution, 

Q(0) = (fi|p|fi), (37) 

(with |f2) being the generalized coherent states for the 
relevant symmetry group) follows from the formal time 
dependence of the density operator 
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H,p 



(38) 



which leads to the following definition of the derivative 
with respect to the dependent parameter: 



d _ 1 
dx\ 2x\ 



^T xV + xV ) s -aif (31) 



To reduce the M independent complex parameters of 
the multimode Heisenberg-Weyl group to the (M — 1) 
independent complex variables parametrizing the SU (AI) 
coherent states, one has to invert the relation between 
the projectors for the multimode Glauber states and the 
SU(M) coherent states \x) N : 



a {a 



E 



N+L i<P(N~L) 

e"l Q l f= |xMx| L 



(32) 



This can be done using the following homomorphism (2lj 



lim ( J~ 

i 2 ^o \ da 



N 



a) {a 



2n 



|x) A r(x| A , (33) 



and the relation 



,2^ + a " 



„2iV 



= Ne~ a 



„2N 



Nl 



(34) 



By means of the relation 







-Q(0,t) =tr(p|fi) (n|), 



(39) 



the properties of the trace and the hermiticity of the 
Hamiltonian one finds 



d_ 

at 



V l (H)-V l (H)*j Q(Sl,t) 
Im (V 1 (H)) Q(n,t), 



(40) 



independent of the specific structure of the dynamical 
group. In the following we will use rescaled units with 

h = i. 

To evaluate the imaginary part of the differential op- 
erator T> 1 (H) for the Bose-Hubbard Hamiltonian (JTJ) , 

" Q{XJ) = -2lJ^e i V\n i )+ l LY.{v\n i f) (41) 



dt 



i=l 



M-l 
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we change once again the parametrization by an ampli- 
tude phase decomposition: 



Xl = VpT, a* = VWe~ i9s 2 < i < M. (42) 



In the case of the Bose-Hubbard model, the pj refer to 
the relative occupation in the j-th well and qj describes 
the relative phase between the j-th and the first well. 

Since the results for the differential operators can be 
used for every Hamiltonian whose dynamical symmetries 
are a SU (M) group, the explicit form of the differential 
operators may be of general interest. A lengthy calcu- 
lation yields the results for the differential operators for 
3 = 1, 



(43) 



k>2 



lm(V l (nif) = Pi ]T 



Pk 



,k.k'>2 



Im (v l {a\a 2 ) +V l {a\a 1 \ 




dp k dq k 



cos q 2 . . ,. 

V Pi ^ 0q k 



, . d 1 [pT d 

+y/PiP2 sin q 2 - h - cosq 2 . — ^— , 

op 2 2 V p 2 aq 2 



and for 2<j<M, 



1 

ImfvUhj)) = --£- 
v v 3 '> 2dqj 



M 



d \ d 

Im(V l (hrf) = Pj \ N -J2p^ )^+Pj 

k>2 



(44) 

d 2 

dp k J dqj ' l ' J dpjdqj 



Im 



d 



(V\a]a 1+1 )+V l {a] +I a, j )) 
= V p-p-T lS m(q j -q j+1 ) (± - — 

1 / \ ( lPi+1 d 

+ 2 C0S W I 



Pi 



<9 



Pi <>'/.; V Pi+1 '''/./ • I 



The advantage of this result is that it directly gives the 
exact phase space dynamics of the Bose-Hubbard model 



in terms of the Husimi- function: 



dQ. 



— (p,q,i) = \ A( + 2^p 2 Pi sinq 2 d P2 



M-l 



+2 ^2 VPk+iPk sin(q k 

Qk+i) (d Pk d Pk+1 ^j 

k=2 

M-l , . 

+ g cosfe +1 - q k )(^ Jj*- 0qk+1 + f-^0 qk )) 

, M M 

+u(n^2( Pi -p k )d qk - ^p k d Pk d qk 

^ k=2 
M 

+ E ( pk ~ Pi)Pk'd Pk ,d, 

k,k'=2 

M d I 

^2( e i- e k)-Q^ | Q(p,q,*), 



fe=2 



k=2 

with the definitions 
qi =0, 



_9_ 



M 

£ 

fe=2 



(45) 



(46) 



Therefore we have derived an explicit formula without 
any approximations. Before we analyze this formula we 
will derive the analogues result for the P-function. 



B. The Glauber— Sudarshan distribution 

The Glauber-Sudarshan or P-distribution is the diag- 
onal representation of the density matrix in the basis of 
the generalized coherent states 

p = j P(Q)|ft) (fi|d/i(ft). (47) 

Since the basis is overcomplete, this description is always 
possible, but not necessarily unique. 

The differential operators for this phase space distribu- 
tion, denoted below by V to avoid confusion, arise from 
a simple integration by parts of the differential operators 
for the Husimi-distribution: 

J P(Q)V l (A)\n) (fi|d/i(ft) (48) 

r v l (A)p(n)\n) (si\dn(Q). 



Thus, one can calculate the time evolution of the P- 
function using the differential operators in an analogous 
way as for the Husimi-distribution: 

P = J P(fi)|ft) (n|d/i(ft) (49) 
= i j (v l (H)* -V l {H)^ P(Q)\Sl) (fi|d/i(fi), 
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or briefly 



d_ 

dt 



P(fl,t) = -21m (v l (H)^j P(fi,t). (50) 



Therefore we can derive the expression for the differ- 
ential operator of the generalized angular momentum op- 
erator by an integration by parts: 



- - d 
V l (E jk ) = -x k —-5 ]k 



(51) 



+x k x*[ (N + M) + -(xV + x*V*) 



In this equation we used the same definitions as above in 
equation (|36[) . The origin of the minor changes compared 
to the case of the Q-function is clear: the additional fac- 
tor M and the 5-symbol result from the different operator 
ordering and the sign is due to the integration by parts. 

Now we can calculate the exact dynamics of the P- 
function for the Bose-Hubbard model with M sites: 



OP, 
~dt 



(P,q,*) = i + A( 2^/p 2 Pi sinq 2 dp 2 



M-l 



+2 ^ VPk+iPk sin(q k - q k +i) (dp k - dpk+i) 
^2 cos(q k+1 -q k )( 1 



k=i 



Pk+i 



dq k +i + J-^uq. 



, M M 

uliN + M^iPi-Pk^qk+Y, 



Pk 

w n2 



Pk 



k=2 



k=2 



dpd 



IV 1 \ 

(Pk-Pi)Pk'd Pk ,d qk J 

k,k'=2 ' 
M 

^(ei - e k )dq k \ P(p,q,t) 

t=9 J 



(52) 



where we used rescaled units h — 1 and the same defini- 
tions (pro| as above. 

A comparison with the result for the Husimi- 
distribution (I45|) shows that due to the operator ordering 
the interaction strength now varies with the particle num- 
ber plus the number of sites, U(N + M). Apart from this 
issue, the first order differential form is exactly the same. 
The second order contribution has apparently the same 
structure as above, but the sign has changed. In both 
cases the second order term vanishes in the macroscopic 
limit N — > oo with UN fixed as 0(1/N). 



C. Liouville dynamics 

In the mean-field limit, the dynamics of a BEC 
in an optical lattice is given by the celebrated dis- 
crete Gross-Pitaevskii equation (GPE) or discrete non- 
linear Schrodinger equation (see, e.g., [251 ] and references 



therein) : 

i±j = ejXj — A(xj+i + Xj-i) + UN\xj\ 2 Xj. (53) 

Using again the decomposition into amplitude and phase 
(|4"2")l . the dynamics can be reformulated as classical 
Hamiltonian equations 



&H 

dpi' 



dqi ' 



(54) 



with the corresponding Hamiltonian function 



M-l 



W(p,q) = -2A Y VPkPk+i cos(<7fc+i - q k ) 



k=l 
M 



M 



-— 2^Pk + 2^ e kPk- 

k=l k=l 



(55) 



One should keep in mind that the parameters of the first 
well are not independent. The GPE describes the exact 
dynamics for vanishing interaction {7 = and an initially 
coherent state, since then an initial state stays coherent 
and the description by a single particle density matrix 
contains no approximations. 

A classical phase space distribution p(p, q, t)dpdq, 
with p, q being canonical conjugate variables, describes 
the probability that an ensemble of particles will be found 
in an infinitesimal phase space element dpdq. The dy- 
namics under the Hamiltonian function TL is governed by 
the classical Liouville equation 



at = dt +{p ' n} = > 



(56) 



where {•, •} denotes the classical Poisson bracket. The re- 
sulting evolution equations for the Hamiltonian function 
(1551 are 



(57) 



dp 



OH dp ^ dTi dp 

dt l~2 9qk dpk k~2 dpk dqk 
= +2A^/p 2 pi sin q 2 d P2 p 
M-l 



+2 A Y VPk+iPk sm(q k - q k +i) (d Pk p- d Pk+1 p) 

k=2 

M-l , \ 

+A Y Qk+1 ~ ■ 



fe=i 

M 



^JPkPk+l 



Pkd qk+1 p + pk+id qk p 



M 



-UNY J {Pi~Pk)d qk p + Y.^- ek ^- d - q - k P 



k=2 



k=2 



A comparison with (|4"5)) and ([52")) shows that the exact 
phase space dynamics consists of a first order differen- 
tial equation plus a many-particle quantum correction of 
second order vanishing in the macroscopic limit N — > oo 
with U N fixed. The first order terms can be thought of 
as the classical evolution since they are identical to the 



8 



results of the Liouville equation. Thus, in the nonintcr- 
acting case this result coincides with the many-particle 
result - the Liouville equation is exact. In this case the 
GPE describes the evolution of the center or maximum 
of the phase space distribution. 

However, the description in quantum phase space goes 
beyond the area of validity of the GPE since there are 
no restrictions on the shape of the initial state, up to 
the usual ones set up by the uncertainty relation. In 
the interacting case the first order part is reproduced 
by the classical Liouville equation, but without the term 
depending on the operator ordering. The first order in- 
teraction term is responsible for a variation of the shape 
of the state, therefore an initial state stays no longer co- 
herent. This fact is usually denoted as the break-down 
of mean-field [26l l27l . [28| , indicating that the description 
by a single mean-field trajectory corresponding to the 
evolution of the center of the coherent state is no longer 
valid. Indeed this breakdown is resolved by using the 
Liouville approach, where we can take into account the 
variation of the shape of the initial state and therefore 
effects due to variation of the higher moments. The sec- 
ond order differential corrections to the classical Liouville 
equation decay with increasing particle number as X/N 
in the macroscopic limit. These terms are responsible 
for many-particle effects as tunneling in quantum phase 
space and (self-)interference. It is interesting to note that 
both the Liouville equation and the whole equation with- 
out approximations conserve the normalization. In a se- 
quel article we will illustrate the methods presented here 
and discuss possible applications [29l |. 



D. Expectation values 

The expectation value of an arbitrary operator B in 
terms of Q- and P-functions is given by the statistical 
average of the phase space distribution 



IB) I P /; (L>)Q(0)<l//(0) 

p(n)Q 6 (n)dfi(n), 



and vice versa: 



(58) 



where P^(f2) and Qg(fl) denotes the (anti)-normally or- 
dered Weyl-symbol of the operator B 

B = J Pg{tt)\n) (fi|d/i(Q) (59) 

Q 6 (n) = (q\b\si). (60) 

In contrast to the symmetrically ordered Wigner func- 
tion the expectation values cannot be expressed in terms 
of one phase space distribution alone. However, the dif- 
ferential algebra formalism allows also to calculate the 
expectation values in terms of the Q-function and the 
differential operator without using the P-representation 



(B) = Tv(Bp) 

= Tr (J B\n) (n\pdfi(n)^j 
v l (B)Q(n)dp(n) 
i> l {B)P(n)dn(n). 



(61) 



Note the interesting correspondence between equation 
(|58[) and equation (jBTjl which reveals the close connection 
between the differential operators and the Weyl-symbols 
of the operator B. 

As an example, we calculate the expectation value of 
the generalized angular momentum operators Ejk — a]afc 
which span the su(M) algebra in the Q-representation: 

jz(Ej k ) = J x k x*Q(x)dp(x) 



~ J £ fe x*(xV+x*V*)Q(x)d M (x) 
y^x*Q(x)d M (x)+0(l). (62) 



and by using the P-function: 

jj{Ejk) = J x k x*P(x)dp(x) 

J -§ J p (x)MV) - - 1 ^P(x)dMx) 

+ 2A7 / ^( xV + x * V *) P ( x ) d M x ) 
= Jx k x*P(x)d»(x)+0(±). (63) 

At the first sight, the expectation values can be decom- 
posed into the classical statistical average and a quan- 
tum many-particle correction that vanishes if the parti- 
cle number N becomes macroscopically large. Moreover, 
we can even concretise the result using an integration by 
parts and the periodic boundary conditions. This pro- 
vides the following result for the Q-function 



(E jk ) = (N + M) j x fc x*Q(x)d Ai (x)-^ fe , (64) 

and the subsequent outcome for the P-function: 

(E jk ) = N J x k x*P(x.)du.(x). (65) 

The differences are of course due to the operator ordering. 
For a coherent state |x ) with P(x) = S(x— x ) we obtain 



(Ejk) = NxkfiX* Q , 



(66) 
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as expected. 



the relative amplitudes and phases 



The calculation of the expectation value of the gen- 
erators of the su(M) algebra by a classical phase space 
average is thus not only a good approximation for large 
particle numbers, but exact. Therefore the only error of 
the expectation values calculated using the Liouville dy- 
namics discussed in section ITV CI is caused by the trun- 
cation of the evolution equations. This error vanishes 
for arbitrary initial states as 0(l/N) in the macroscopic 
limit TV — > oo with UN fixed. 



THERMODYNAMICS 



The method presented above is not restricted to an 
analysis of the time dependence of the system, there are 
multifarious applications. As an example, we consider 
the (unnormalized) density operator of the canonical en- 
semble, 



-I3H 



(67) 



with (3 = 1/kT and h = 1. This expression describes 
the quantum mechanical version of the canonical parti- 
tion function in statistical mechanics obeying the Bloch 
equation 



dp 
dp 



■(pH + Hp). 



(68) 



A. Thermodynamics of the Q— function 



<9Q(p,q) 



dp 



M 



k=l 



M 



^2 e kPkPk'd Pk , - ^2 e k Pkd Pk 

k=l,k'=2 k=1 
M-l , 

+A ^ ( 2^/pkPk+i cos^+i - q k ) 

k=l ^ 

M v 

x ( N -J2 Pkldpk/+ d Ph+1 )\ 

k'=2 ' 

M-l i 

+ T \ — — sm(q k+1 -q k )(d qk - d qh+1 ) 
2 ^ V Pk+i 

M-l 

+A ^ VPkPk+i cos(q fe+ i - qk)d Pk 



k=2 



UN(N - 1) 



M 



M M 



k=l 



k=l k'=2 



M 



M 



uY,pU(N-i)-J2p»' d p k ') d p k 



u 



k=2 k'=2 
M M 

J2 P k Pk'Pk"d Pk ,d Pk „ 
k=l k'.k"=2 



M 2 



k=2 
M 



k=2 



■g £ d qk d qk ,\Q(p,q) 

k,k'=2 



(70) 



Besides the lengthy expression one already recognizes an 
underlying structure: The leading terms of each contri- 
bution show a close analogy to the GPE Hamiltonian 
function Before we have a closer look at the con- 

nection to the classical result, we derive an expression 
for the solution of the Bloch equation in terms of the 
Glauber-Sudarshan distribution. 



Translating the relation (|68j) into differential operators 
acting on phase space densities, namely the Q-function, 
yields the formal result: 



dQ 

dp 



-Re 



(V\H)) Q. 



(69) 



B. Thermodynamics of the P— function 



Analogously to the case of the Husimi-distribution one 
can derive the result for the P-function, 



Analogously to the calculations in the previous section, 
we can evaluate the real part in the parametrization of 



dP 
~dp 



-Re 



(VI{H) 



(71) 
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Here the evaluation of the real part yields 



d-P(p,q) 

8(3 



A I 



A I 



= \ -(N + M)J2 ( -kPk + ek 



At 



k=l k=l 
A I 



k=l 



^2 e kPkPk'd Pk , + 2J £kPkd Pk 

fe=l,fe'=2 k=2 
M-l , 

A ^ ( ^PkPk+l COs(q k+1 - q k ) 
M 

((N + M)-^2Pk'd Pk ,+-d Pk+1 

y 



k>=2 



M-l 



A \ - / Pk 

' 2 ^ V p fe+ i 



k=l 
M-l 



sin(g fe+ i -q k )(d qk - d qh+1 ) 



-A 2J v'PkPk+i cos(g fc +i - 



fc=2 



,1/ 



^ + M)(iV + M +1 ) I:p| + a(2Ar + M) 

fc=l 

M 



fc=2 



+2C/ ^ Pkpyd Pk , - ^J2 p k E Pk>Pk»d Pk ,d Pk „ 



k=l,k'=2 
AI M 



fc=l,fc'=2 
M Af 



fc=l k',k"=2 



k=2 k'=2 
M 



?E^ + ?EC 



fe=2 
At 



k=2 



+ g E P(p,q). 

fc,fc'=2 J 



(72) 



Note the subtle, however important differences due to 
operator ordering compared with equation (|70[) . 



C. Classical vs. quantum statistical mechanics 



The distribution function of the classical canonical en- 
semble 



p = e 



-PNH 



(73) 



given in terms of the Hamiltonian function (|55[) solves 
the Bloch equation 



dp 
dp 



-Hp 



{IVl — J. 
2AiV ^2 VPkPk+i cos(q k+1 - q k ) 
k=l 

TIN 2 M M 1 

-^Eri-* ( 74 ) 
fc=i fe=i j 

A comparison with equation (|7D|) and equation ([72")) 
shows that the quantum many particle Bloch equation 
and its formulation in terms of the Q— and P-function 
can also be separated into a classical contribution, the 
leading order of N, which is governed by the Gross- 
Pitaevskii Hamiltonian function and quantum correc- 
tions. Amongst others these additional quantum terms 
ensure the minimal uncertainty for low temperatures. 
The high temperature limit is in both cases given by an 
uniform distribution. 



VI. CONCLUSION AND OUTLOOK 

In this paper we have developed phase space techniques 
which provide an alternative tool to investigate and ana- 
lyze the dynamics of one-dimensional M-site, TV-particle 
Bose-Hubbard systems. The quantum phase space is con- 
structed in terms of generalized SU (M) coherent states, 
which conserve the number of particles. This changes the 
corresponding phase space to a compact manifold. In the 
context of Bose-Einstein condensates, the SU{M) coher- 
ent states have a special significance for these systems as 
they describe fully condensed states. 

The phase space dynamics can be treated efficiently in 
terms of the differential algebra developed by Gilmorc. 
In this way the su(M) operator algebra is mapped onto 
differential operators acting on the multimodc coherent 
states. The resulting evolution equations for the (gen- 
eralized) Husimi (Q) and Glauber-Sudarshan (P) phase 
space distributions are second order differential equa- 
tions. These (exact) evolution equations provide a con- 
venient starting point for further developments. 

Firstly, it is immediately observed that the second or- 
der terms scale as 1/N and therefore vanish in the macro- 
scopic limit N — > oo with UN fixed. For large N, the 
evolution reduces to first order equations of the form of 
(classical) Liouvillian dynamics. The phase space ap- 
proach therefore provides a remarkable direct derivation 
of the celebrated many-particle mean-field limit. 

Secondly, this phase space method offers a clue to 
generalize the mean-field approximation, which describes 
strongly localized quantum states by a single point in 
phase space. Arbitrary quantum states can be repre- 
sented by an ensemble of phase space trajectories, which 
is constructed to approximate the initial quantum phase 
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space (Husimi) distribution. Then each trajectory fol- 
lows the (classical) mean-field equations. This allows a 
straightforward computation of expectation values. 

Thirdly, the resulting second-order partial differential 
can be attacked directly by numerical methods. 

Finally, there is the challenge to explore the regime 
between the (classical) mean-field description and the 
full quantum dynamics by generalizing the semiclassical 
phase space methods developed during the last decades 
for the flat space to systems with SU(M) symmetries and 
a compact phase space. 

In addition, it should also be noted that the evolution 
equations can also be generalized to master equations 
describing systems coupled to an environment or systems 



with an effective decay. In future work we will address 
some of these problems, starting with first applications 
to the two-mode Bose-Hubbard system in a forthcoming 
article [H. 
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